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Abstract 

Using the graded eigenvalue method and a recently computed extension of the 
Itzykson-Zuber integral to complex matrices, we compute the fc-point spectral cor- 
relation functions of the Dirac operator in a chiral random matrix model with a 
deterministic diagonal matrix added. We obtain results both on the scale of the 
mean level spacing and on the microscopic scale. We find that the microscopic spec- 
tral correlations have the same functional form as at zero temperature, provided 
that the microscopic variables are rescaled by the temperature-dependent chiral 
condensate. 
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1 Introduction 



In the past few years, random matrix models have started to be used in or- 
der to understand details of the spectral properties of the Dirac operator of 
non-abelian gauge theories Since random matrix models are free of 

dynamical input, one expects them to describe universal spectral properties, 
i.e., properties which do not depend on the specific details of the dynamics. 
In the case of non-relativistic quantum mechanics, it has long been known for 
a variety of physical systems that random matrix models are able to describe 



the universal spectral correlations on the scale of the mean level spacing p2 



p4| . In the case of the Dirac operator of QCD the chiral structure has to be 
incorporated, requiring the introduction of the chiral random matrix ensem- 
bles |Q . The symmetry constraints lead to the definition of the chiral Gaussian 
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Orthogonal (chGOE), Unitary (chGUE) and Symplectic (chGSE) Ensembles. 
A precise classification can be found in Ref. |7|. Our chiral random matrix 
model is based on unitary symmetry and will be defined in Sec. |2|. 

For the case of QCD with two colors, the universality of the spectral corre- 
lations in the bulk of the spectrum of the lattice Dirac operator was demon- 
strated convincingly in Ref. where predictions from random matrix theory 
were compared with lattice data obtained by Kalkreuter. However, the fact 
that the chiral symmetry of the Dirac operator is spontaneously broken im- 
plies that another type of universality exists concerning the microscopic limit 
of the spectral density. The spectral density is defined by 

p(\) = (£5(\-\ n )) A , (1) 

n 

where the \ n are the eigenvalues of the Dirac operator iD = i^^d^ + J^A^ 
(in euclidean space) and the average is over all gauge field configurations A 
weighted by the euclidean QCD action. Note that iD is hermitian so that the 
\ n are real. Moreover, because iD anticommutes with 75 all eigenvalues occur 
in pairs ±A so that p is an even function. The chiral condensate H is related 



to the spectral density via the Banks-Casher formula [25 



H = £p(0) , (2) 
where V is the space-time volume. The microscopic spectral density is defined 

by § 

1 u 

Ps(u) = lim — p{— ) . (3) 

V — >oo y c V c 



Based on sum rules derived by Leutwyler and Smilga p6| , p s was conjec- 



tured to be a universal quantity |§. This conjecture was supported by several 
works P,P, p!3| , p!9| , |27| -p9| . As a universal function, the microscopic spectral den- 
sity can be computed in the framework of random matrix models. This task 
was performed at zero temperature for all three chiral ensembles [|J-|[J. The 
most direct evidence for the universality of p s came from a recent comparison 
of lattice data obtained by Berbenni-Bitsch and Meyer with predictions from 
random matrix theory where impressive agreement was obtained ||30||. 



Recently, the chiral random matrix models were extended to finite temperature 



T and finite chemical potential p ^DQ-T$\. This leads to the addition of a 



deterministic diagonal matrix to the random matrix used at zero T and p. 
In such models, one can investigate, e.g., how the chiral condensate decreases 
as a function of T and/or p. For the simplest model in which the additional 
diagonal matrix is given by the lowest Matsubara frequency ttT and, thus, 
proportional to the unit matrix, the spectral density has been computed in 
Refs. p3|,|T9! and its microscopic limit in Ref. [[Til. In the present paper, we 
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(i) generalize these results to arbitrary diagonal matrices and (ii) compute 
all higher spectral correlation functions, both on the scale of the mean level 
spacing and in the microscopic limit. The first extension is important since 
it is often desirable to consider more general diagonal matrices, e.g., if one 
wants to include higher Matsubara frequencies or in the approach of Ref. Jl"5|| 
where the formation of instanton-anti-instanton molecules was modeled. The 
second extension allows for further analysis of lattice results for, e.g., scalar 
susceptibilities. Our results also add further evidence to the conjecture that 
the microscopic spectral correlations are universal. 

This paper is organized as follows. In Sec. ^|we define the model and summarize 
some known results which will be refered to in the following. In Sec. |3| we set up 
the computation of the spectral correlation functions in the framework of the 
graded eigenvalue method and express the result in terms of a double integral. 
This double integral is evaluated in Sec. [| for the spectral density and the 
spectral correlations on the scale of the mean level spacing and in Sec. [5| for the 
microscopic limit of the spectral correlation functions, respectively. The results 
are discussed in Sec. |6|. Three appendices are provided to discuss convergence 
issues and technical details associated with saddle-point approximations in 
Sec. |. 



2 Formulation of the problem 

We formulate the problem in the framework of random matrix theory in 



Sec. |2.1| . In Sec. [2.2| , we discuss the connection with Dyson's Brownian motion 
model. 



2. 1 Definition of the model 



We are interested in studying the spectral correlations of the Dirac operator 
on the scale of the mean level spacing and in the microscopic limit. Since the 
spectral correlations are universal in these two limits they can be obtained in 
the framework of a random matrix model ||. In such a model, the matrix of the 
Dirac operator (in euclidean space) is replaced by a random matrix with ap- 
propriate symmetries, the average over gauge field configurations is replaced 
by the average over random matrices, and the complicated weight function 
containing the euclidean action is replaced by a simple Gaussian distribution 
of the random matrices. Moreover, at finite temperature a deterministic di- 
agonal matrix with real entries has to be added to the random matrix of the 



Dirac operator |l^-[l^,p^| . (Note that one could also add a non-diagonal de- 



terministic matrix since it could always be diagonalized by a suitable rotation 
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of the basis states. This is possible due to the invariance of the random matrix 
ensembles with respect to basis changes.) 



The matrix representing the Dirac operator in a chiral basis has the structure 



A 



imll N W + Y 
W ] + Y imH- N 



(4) 



where m is an additional mass term, W is a square random matrix of di- 
mension N, and Y = diag(yi, . . . ,un) is a diagonal matrix with real entries. 
The total number of eigenvalues of A is 2N which can be identified with the 
space-time volume V in Eqs. (|2|) and (|3|). The specific form of Y is model de- 
pendent, but we will obtain results for an arbitrary choice of Y. Note that we 
consider the quenched case where Nf, the number of flavors, is equal to zero. 
Furthermore, we consider only square matrices W. Using rectangular matrices 
one can model the topological charge v which corresponds to the difference 
between the number of rows and columns of W. The extension of the present 
calculation to non-zero Nf and v is in progress. 

In this paper, we study the chiral GUE appropriate for QCD with three or 
more colors for which W is a complex matrix. The probability distribution of 
W is a Gaussian, 

P(W) = f — — I exp (-NY?tiWW ] ) , (5) 

where S is a real parameter equal to the chiral condensate at zero temperature. 

We shall consider separately the case where all the y n are equal, i.e., where 

Y = yUjv (this will be refered to as the "special case") and the case where 

Y is arbitrary (this is the "general case"). This separation is useful since the 
various mathematical structures are most easily identified in the special case 
from which they can then be generalized. 

The chiral condensate 2, which depends on the matrix Y, has already been 
computed for both the special case [[14]], 



2 = EJl - (m 2 , (6) 



and the general case jl5i , 

2 = Y?x . (7) 
Here, x is the only real and positive solution of 
2 1 N 1 

" 2 = nE^T¥ (8) 
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or zero if no such solution exists. (See Appendix y for a proof that this equation 
has at most one real and positive solution.) 

The eigenvalue density and its microscopic limit for the special case were 
computed in Ref . |L9| . In the following, we will compute the eigenvalue density 
and the spectral correlation functions on the scale of the mean level spacing 
as well as the corresponding microscopic limits for the general case. 

The /c-point correlation functions measure the probability of finding ener- 
gies around the points x\,...,Xk, regardless of labeling. Apart from trivial 
5-functions they are defined by 

Rk{xi, ...,x k )= (--)" I d[W]P{W) fl Im tr^— , (9) 

where Xp = x p + ie with e positive infinitesimal. In particular, we have p(X) = 
i?i(A). It is technically simpler to compute the functions 

Rk{xi, ...,x k )= (--)" f d[W]P{W) fl ^^Z^ (10) 

\ 71/ J p=1 X p A 

from which the R k can be reconstructed by taking the imaginary parts of the 
traces. 



2.2 Joint probability density function of the radial coordinates 



Dyson's Brownian motion model of random matrix theory [^[H] describes 



the transition from an arbitrary to a Gaussian ensemble. Since this approach 
is often used, we establish the link to our random matrix model. The main 
modification is of course due to the chiral structure. In Ref. J^j, we gave a 
detailed discussion of a diffusion equation in the space of arbitrary complex 
ordinary and supermatrices. In the case of ordinary matrices, this diffusion is 
the chiral analogue for the chGUE of Dyson' Brownian motion for the GUE. 

The pseudo-diagonalization of the matrix W can be written as W = UXV, 
where U e U(N), V G U(N)/XJ N (l), and X = diag(xi, . . .,x N ). The x n are 
referred to as radial coordinates, they are defined on the positive real axis. 
The transformation of the Cartesian volume element to radial and angular 



coordinates reads 32 



N 

d[W] = A 2 N (X 2 ) J] x n dx n dp(U)dp(V) , (11) 

n=l 

where Ajv(X 2 ) = Y[ ri<m i.x 2 l — x^) is the Vandermonde determinant and dp{U) 
and dpiV) are the invariant Haar measures. 
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Consider now Y = in Eqs. (f|) and Obviously, the integrand depends 
only on the radial coordinates and the group integrals are trivial. This so- 
called rotation invariance is broken for non-zero Y. Nevertheless, even in this 
case, the mathematical problem can be reduced to integrals over the radial 
coordinates. By shifting the real part of the diagonal of W and thus also of 
W> by the matrix Y, the breaking of the rotation invariance is removed from 
the traces in Eq. (Q). The resulting probability density function P(W — Y), 
however, now depends explicitly on the angular coordinates U and V. To 
reduce the problem to the radial coordinates, one has to perform the group 
integrations 

T(X,Y) = J dn(U) J d/i(V)P(UXV - Y) (12) 

which can be viewed as the extension of the Itzykson-Zuber integral to complex 
matrices. In Refs. [32^D|, this integral was calculated, 

-rrv v \ _ 1 det[7(x n ,y m )] 

n,m=l,...,N /, q\ 

UA ' yj = M A N (X*)A N (Y*) ' [L6) 
where the entries of the determinant are given by the functions 

7 (x n , y m ) = 2NT 2 exp (~NY?(x 2 n + y 2 j) I (2NE 2 x n y m ) . (14) 
Here, Io(z) is the modified Bessel function of zeroth order. 

Collecting everything, the /c-level correlation function is obtained by integrat- 
ing the distribution 

N i a (X 2 ) N 

T(X,Y)A 2 N (X 2 ) J] x n = det[ 7 (x n ,y m )] n , m=1 ,...^-^-( J[ x n (15) 



over the N — k radial coordinates x^+i, . . . , x^. The expression ( |15|) is called 
the joint probability density function. Thus, we have identified the kernel of 
the diffusion equation of Ref. with the joint probability density function 
of our random matrix model. Apparently, it is highly non-trivial to actually 
perform the integrations over the N — k radial coordinates. In Sec. |3.2| , we 
will show that the super symmetric analogue of the joint probability density 
function can be used for a much faster and explicit calculation of the spectral 
correlation functions. 



3 Correlation functions and supersymmetry 



After the introduction of random matrix theory by Wigner [33], many of the 



mathematical details were worked out by Dyson pTf] and Mehta [53]. In 1983 
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Efetov introduced the supersymmetry method |33] which was further devel- 
oped in Ref. |35 |. In Sec. |3j] we employ these methods to construct a supersym- 
metric representation for the spectral correlation functions. It turns out that 
the formal aspects of the problem discussed here are related to the transition 
from Poisson regularity to the GUE. In Ref. |JB), this transition was discussed 
with supersymmetry by extending the graded eigenvalue method developed in 
Ref. [37|. In Sec. |3.2| we further generalize this approach to the chGUE. The 
main ingredient is now the supersymmetric extension of the Itzykson-Zuber 
like integral pil mentioned in the previous section. 



3.1 Supersymmetric representation 



The correlation functions 
function 



and (|ToD can be obtained from the generating 



f r , / x n^-i det(x+ + J v - A) 



(16) 



where J stands for J 1; 
by 



R k (xi, 



Jk- Note that Zk(0) = 1. The R k are then generated 



X k) = -7T- 



1 

2~n 



■z k {J) 



(17) 



The ratio of the products of determinants in Eq. flIE| ) can be written as an 
inverse graded determinant which, in turn, can be expressed as an integral 
over graded vectors. We define ^-dimensional vectors z pl , z p2 , Cpi, (p2, and 



2A^-dimensional vectors z p = (z p i,z P 2) and ( p = (Cpi, Cp2) 



The z p contain 



commuting variables whereas the ( p contain anticommuting variables. These 
variables are combined in a graded vector ip = (zi, . . . , Zk, Ci> • • • > Ck) T of di- 
mension 4kN. The generating function (|T6|) can then be written as 



Z k (J) = J d[W}P(W) detg^D 

d[W]P(W) f d[tfj} exp (i^Dtfj) 



18) 



where D = (x + + J) <S> 2-2AT — ^-ik ® A is a matrix of dimension AkN. Here, 
x + J = diag(xi — Ji, . . . , x k — Jk, %i + Ji, ■ ■ ■ , x k + Jk) is a graded matrix 
of dimension 2k, and our convention for the direct product of two matrices 
S and T with number of rows r$, and number of columns c$, Ct is that 
(S (g> T)(i_iy T+ k,{j-i)c T +i = SijT k [. The integration over W can now be per- 
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formed. We write 
with 



W 
t0 



V> = tr (WE + W^tf 



(19) 



f — y^(^p2^ D i — Cp2c p i) 
P =i 

and obtain the ensemble average 



(20) 



/ d[W] P(W) exp (-Ht(WE + W ] E ] ) 



exp 



1 



-tiEE ] 



exp 



1 



rtrgF!F 2 . 



NY? J * \ NT? 

The Fi {i = 1, 2) are 2/c-dimensional graded matrices of the structure 

pBB pBF 

i i 
pFB pFF 



(21) 



[22) 



where the boson-boson and fermion-fermion block both have dimension k. The 
entries are given by scalar products 



( F BB\ = t 
(■^i )pp' — Z piCp'i 



(■^i )pp' Cpi z p'i 
(-^i )pp' ~~ CpiCp'i 



(23) 



with = 1, . . . , fc. We now perform a Hubbard-Stratonovitch transforma- 
tion to decouple Fj and F 2 , 



exp 



rtrgFiFa 



2 2fe2 J d[a] exp (-iVE 2 trg(7(7 t + i trg(F l( x + 



(24) 



where a is a complex graded matrix of the same structure as F^. Combining 
Eqs. (HD, (0), and (H), we obtain 



Z fe ( J) = 2 2k2 J d[ip] J d[a] exp (-NY?tigoa ] + z'^Gty 

with 

„ H /v + 




Uiv 



+ (x + +j) ® u 2 jv - n 



2fe 



o y 

Y 



(25) 



.(26) 
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Note that we have set the quark mass m equal to zero in (^6]). It can always 
be reinstated at the end of the calculation if we shift the x p by im. 

We would now like to interchange the order of the a- and the -^-integrations 
in (p5|) and to perform the integration over if). This would be allowed if the if>- 
integral were uniformly convergent. Uniform convergence is guaranteed for the 
ordinary GUE but not for the chGUE. In fact, the ^-integral as it stands is not 
uniformly convergent. This presents certain problems which are discussed in 
some detail in Appendix There are various ways to proceed. For the present 
calculation, we note that there is a very strong conjecture that ignoring the 
convergence problem will lead to a result whose real part is incorrect but 
whose imaginary part is correct. Since we are only interested in the spectral 
correlation functions which are given by the imaginary part, we choose to go 
ahead with the calculation without modifications. However, the reader should 
consult Appendix [A| for a more complete discussion of this problem. 

Ignoring the convergence issue, we now shift a and ert by — (x + J) and perform 
the integration over if) to obtain after some algebra 

Z k {J) = 2 2k " J d[a] exp (-iVS 2 trg( ( r - (x + J))(cr t - (x + J))) 

x detg- 1 (a + (a + y <g> ti N - l 2fc ® Y 2 ) . (27) 



In order to work out the integrals over the supermatrix, there are various 
ways to proceed. In particular, Efetov's MM coset method has been used in 



numerous applications. It is based on a saddle-point approximation in the 
large- N limit. To solve the problems outlined in the introduction we use the 
graded eigenvalue method |]37|j36| which is, in this context, particularly well 
suited. 



3.2 Reduction to eigenvalue integrals 



The matrix a can be written in pseudo-diagonal form as a = usv, where u G 



U(fc/fc), v G V(k/k)/V 2k (l), and s = diag(si,is 2 ) 



with Sj 



diag(sij, 



Skj) 



for j = 1, 2 [32j| . The s P j are real and non-negative. The transformation of the 



Cartesian volume element to radial and angular coordinates reads [32 



d[a] = J(s)d[s]dfi(u)dfi(v) 



d[s] = Y[ ds p ids P 2 
P =i 

J(s)=B 2 (s 2 )Y[ Spl s p2 
P =i 



9 



with 



B k (s) 



AJx) 



A k (si)A k (is 2 ) 
U P , q (s pl - is q2 ) 

■ JJ (x p — Xq) . 

p<q 



The group measure is the invariant Haar measure. We obtain 



(29) 
(30) 



Z k (J) = 2 2fc2 exp (-iVS 2 trg(x + If 

x j d[s]J{s) exp(-iVS 2 trgs 2 )detg^ 1 ((s + ) 2 <g> t N - t 2k ® Y 2 

dfx(u) J dfx{v) exp (iV£ 2 trg(a + a r )(x + J)) . (31) 



The integral over the diagonalizing groups has been computed in Ref. |32 
The result is 



J dfi(u) J d/j,(v) exp (iV£ 2 trg(<7 + a^)(x + J) 

_ {2NT 2 ) 2k det[I (2NZ 2 s pl (x p ,- J p ,))} det[J (2iVS\ 2 (3 V + J q ,))] 



2 2k \k\y 



Bk(s 2 )B k ((x + J) 2 ) 



(32) 



with p,p' = 1, . . . , k and q,q' = 1, . . . , k. Here, Jo and Jo are the Bessel and 
modified Bessel functions of order 0, respectively. Thus, Eq. (pTf) becomes 



„ , , , (2iVS 2 ) 2fc exp(-iVS 2 trg(a; + J) 2 ) , ,. . „ 
Z k (J) = 1 + - /m9r> 77Z , T \Z\ — / d l s ] II s pi s p2B k (s 2 ) 



P =i 



(k\yB k (( x +jy) 

x exp(-iVS 2 trgs 2 )detg" 1 ((s+) 2 ® 1 N - TL 2k <g> F 2 ) 
x det[/ (2A^E 2 Spi(x p / - Jp'))W=i,-,& 
x det[J (2A^S 2 s 92 (a; (? ' + J g '))W=i,..,fc • 



(33) 



Note the contribution of unity to the generating function. It ensures the correct 
normalization of Z k at J = and is due to so-called Efetov-Wegner terms 
which were discussed in detail in Ref. |32 . 



The calculation now greatly simplifies due to the determinant structure of the 
function B k (s). In fact, Eq. (|2D|) can be rewritten as 



B,(s) 



det 



Spi — is q2 



(34) 



p,q=l,...,k 
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Performing the differentiations with respect to the J p according to Eq. (pj]) and 
using the determinant structure of Z k , straightforward manipulations show 
that the /c-point functions can be written as 

R k {xi, ■ ■ ■ , x k ) = det[C N (x p , ^ 9 )]p, ? =i,... )fe , (35) 

where the function CV is given by 



C N {x pi x q ) = -{2NT?) 2 x p I I d Sl ds 2 ^ 2 exp (-iVS 2 (s? + si)) 

IT J J Si ~\~ Sn 

Li 

N v 2 + s 2 

xI (2NE 2 x pSl )J (2NE 2 x q s 2 ) Im ^7 • (36) 

Hence, all fc-point functions are known if the double integral ( |5ED can be 
computed. This can be done by saddle-point approximation in the large- N 
limit. However, for some applications it is more convenient to rewrite the 
double integral as an integral over a 2 x 2 complex graded matrix a by using 
our result ( |32"1) for the integral over the diagonalizing groups in the opposite 
direction. This is by no means a trivial step since we needed the angular 
integral in order to obtain the determinant structure of fl35[) . Using Eq. (|32|) 
for k — 1, we obtain from ([36]) 



8 e N ^ x p- x ^ f ( t x 

C N (x p ,x q ) = -x p — 2 - 2 Im / d[a] exp (-N£(a, a T )J (37) 

C(a, a ] ) = S'trgaa 1 " - ^trg log ((a + + s)((a+) t + x) ®H N -H 2 ®Y 2 ) , (38) 

where a is now a 2 x 2 complex graded matrix and x = diag(x p , x q ). Note that 
a and have been shifted by x to obtain Eqs. ([37D and d38|) . The integral in 
Eq. (|3"7D can also be computed in saddle-point approximation in the large- N 
limit. Since the matrix a can be pseudo-diagonalized, a = usv, we can assume 
that a is diagonal at the saddle point, a = diag(si, is 2 ), and thus obtain from 
(I31 



where we have omitted the imaginary increments. 

We will use Eq. (|37| ) to compute the spectral density and the universal spectral 
correlations in the bulk of the spectrum in Sec. || and Eq. (^) to compute 
the microscopic limits of the spectral density and the spectral correlation 
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functions in Sec. [)} Our expressions allow for a relatively straightforward finite- 
N analysis. However, we are eventually interested in the thermodynamic limit. 
Therefore, we compute our results in the limit iV — > oo. 



4 Spectral density and universal spectral correlations in the bulk 
of the spectrum 



In this section, we compute the spectral density and the universal spectral 
correlations in the bulk of the spectrum on the scale of the mean level spacing, 
for both the special and the general case. The mean level spacing D(x) is 
defined as D(x) = l/p(x). With our definitions, we have D(x) ~ 1/(JVS) in 
the bulk of the spectrum. Thus, the correlations are to be computed for energy 
differences ~ l/(iV£). We discuss the special case in Sec. fO| and turn to the 
general case in Sec. [O]. 



4-1 Special case 



In the special case Y = yH N with y real, (|39|) becomes 

C( Sl , lS2 ) = E 2 ( S 2 + si) - log { j S2 + Xq £~ y2 2 ■ (40) 

(si + x P Y - y 2 

To minimize C, we differentiate, first with respect to si, 

^ = 2 E % + r ^i±^i0 (41) 
dsi (si + x p y - y 2 

which yields a cubic equation for si. Differentiation of C with respect to is 2 
yields an identical equation with s\ and x p replaced by is 2 and x q , respectively. 
To compute the spectral density, we have to take the limit To compute 

the universal spectral correlations, Xp Xq IS of order l/(iVE) and, thus, x p = x q 
in the large- N limit. Hence, the saddle-point equations yield identical solutions 
for s\ and is 2 , and we define si = is 2 = s at the saddle point. Furthermore, we 
define x = (x p + x q )/2 and Ax = x p — x q . Note that x p = x q = x for N — >• oo. 
The saddle-point equation then reads 

s 3 + 2s 2 x + s(l/S 2 + x 2 - y 2 ) + xjY? = . (42) 

The Lagrangian ( fHf ) at the saddle point can then be expanded in Ax, 

, (s + x - Ax/2) 2 - y 2 2(s + x) . o 

Co = - log ) l} \ w K ' 2 Ax = -2E 2 sAx , 43 

(s + x + Ax/2) 2 — y 2 [s + x) — y 2 
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where we have used Eq. fl4~l|) in the last step. The integral over the quadratic 
fluctuations yields a factor of 1/4 (cf. our convention for the normalization of 
the Gaussian in Eq. (p4[)), and we obtain from Eq. ( p7|) 



O O 1 H -L 

C N (x p , x q ) = -x p — — - Im exp (-NC ) 



-x\ 4 

Im exp(2iVE 2 sAx) 



71 X^ — X 2 q 4 

e 2NY?xAx 



71 Ax 

exp(2NE 2 (x + Re s) Ax) . (44) 



sin(2iVS 2 Aa;Ims) ,, 1V ,. 2 



ttAx 

Because of the identity 

det [f{x p , x q ) exp {a{x p - x q ))] p q=1 ^ k = det [f{x p , x q )} p q=1 ^ k (45) 
for arbitrary / and a we obtain in the large- N limit 

R k (x) = det[K N (x p ,x q )] M=1> ... tk (46) 

with 

sin(2iVS 2 (x p -x g )Img) 

K iV (Xp,XgJ = — r . (47) 

71 [Xp Xq) 

The spectral density now follows immediately, 

2iVS 2 

p(x) = K N (x, x) = Ims. (48) 

71 

The cubic equation ([42"D always has one real solution. If the other two solutions 
are real, K N (x p ,x q ) and p(x) are zero for this value of x. Otherwise, the 
other two solutions are complex conjugates of each other, and we have to take 
the solution with positive imaginary part because of the positive imaginary 
increment of s. This is the solution of ( f4"2"D which enters ( |47| ) and ([48]). 

Furthermore, we obtain from ([48"1) 

D(x) = Svfbi (49) 

so that ( pETf ) can be rewritten as 

^^^ s 't^" (50) 

Eq. (^) reproduces the well-known result for the spectral correlations on 
the scale of the mean level spacing. Thus, neither the chiral structure of our 
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random matrix model nor the addition of the diagonal matrix Y changes these 
universal correlations. 



We would like to emphasize that the result for p(A) has been obtained previ- 
ously in a simpler way |16| , 19 . 18| . However, our method has the virtue of also 
yielding all higher-order correlations and their microscopic limits in one single 
formalism. 



4-2 General case 



Consider now the case of general (but real) Y where the Lagrangian has the 
form (|3D|). Again, we minimize £ and obtain identical solutions for si and 
is2 in the limit x p = x q . With the definitions of the previous subsection, the 
equation for s is now 

1 N s + r 

This is a polynomial equation of order 2N + 1. We now demonstrate that 
this equation always has at least 2N — 1 real solutions. The remaining two 
solutions are either real or complex conjugates of each other. 

Define z = s + x and consider the function 

1 N z 

/(,) = S 2 (, -*) + -£-—• (52) 

iv n =l 41 »n 

First assume that all the y n are different and non-zero. The function f(z) has 
2N poles at ±y n (n = 1, . . . , N). Now consider two adjacent poles z m < z m+ i. 
We have f(z m + e) > and f(z m+ i — e) < 0, where e is positive infinitesimal. 
Since f(z) is continuous for z m < z < z m+ i, there must be a zero of / in this 
interval. There are 2N — 1 such intervals, hence there are 2N — 1 real zeros of /. 
This proves our assertion. If some of the y n are equal or zero, the saddle-point 
equation (|5T|) is of correspondingly lower order, and our argument applies as 
well. 

In analogy to the result of the previous subsection, only the (single) saddle 
point with positive imaginary part contributes to CN(x p ,x q ). Expanding the 
Lagrangian about this saddle point yields 



_ 1 f Js + x-Ax/2y-yl ^ 2 » (g + x) 
g (s + x + Ax/2y-yl~ N ^ (s + x) - y* 

= -2T?sAx (53) 
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using ([51]) in the last step. Going through the same steps as in the previous 
subsection, we obtain the same results as before, the only difference being that 
s is now the only solution with positive imaginary part (or zero) of the more 
complicated equation fl5lD instead of the cubic equation (|42|) . 



5 Microscopic limit 



In this section, we compute the /c-point functions in the microscopic limit, i.e., 
for x p ~ l/(iV£). In Sec. |Q| , we separate the radial integrals in (|36|), then 



consider the special case Y = yTL N in Sec. |57|, and finally turn to the general 
case in Sec. [Oj Although the special case obviously follows from the general 
case, we find it useful to discuss the special case first in order to make the 
structure of the radial integrals more explicit. The final result will be that 
the fc-point functions are given by the zero-temperature result provided that 
the microscopic variables are rescaled by the temperature-dependent chiral 
condensate. 



5. 1 Separation of the radial integrals 



We first define microscopic variables u p = 2Nx p and u q = 2Nx q . The two 
integrals in Eq. (^) can be decoupled by writing 



exp(-NZ 2 (s 2 + s 2 )) 
s\ + sl 

We thus obtain 



NY? J dt exp {-NT?t(s\ + s 2 2 )) . (54) 



4iV 2 £ 6 x 



C N (u p ,u q ) = u p / dtf 1 (t)f 2 (t) (55) 

7T J 
1 

with 

00 N 1 

f 1 (t)= [ ds lSl exp (-NZHsj) lo^VO Im ft 2 _ 7 ~ ry 2 (56) 

^ n=l Vn l s l + le ) 

oo 

Mt) = J ds 2 s 2 ex V (-NZHsl) J (Z 2 u q s 2 ) f[{y 2 n + s 2 2 ) . (57) 
o n=1 

These two integrals can now be computed in saddle-point approximation. 
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5.2 Special case 



Consider the special case Y = yH-N with y real. The critical value of y above 
which the chiral condensate vanishes is given by y c = 1/S so that we restrict 



ourselves to the interval < y < y c in the following. Eqs. ( p6|) and fl57|) become 



oo 

/i(*) = / ex P (-NZHsj) I (E VO Im ^ 2 _ ^ - . £)2)JV (58) 

oo 

h{t) = J ds 2 s 2 exp (-NZHsl) J Q (Y?u q s 2 ){y 2 + s 2 2 ) N . (59) 
o 

These two integrals can be evaluated in saddle-point approximation in the 
large- iV limit. This somewhat technical task is performed in Appendix [FJ. It 
turns out that depending on the value of t, the product of f\(t) and f 2 (t) can 
be of different order in 1/N. Specifically, there is a "critical" value of t, 

«• = <m> • (60) 

below which fi(t)f 2 (t) is of order 1/N and above which fi(t)f2(t) is suppressed 
exponentially. Note that t c > 1 since y < y c . Thus, only the interval 1 < t < t c 
contributes to the integral in Eq. ( |55|) in leading order in 1/N. We obtain from 
Eqs. O and ( pT5|) 



7T 



/i(*)/ 2 (t) = Jb ( EupVV* - V*c ) Jb ( Eu fl VV* - V*c ) (61) 



to leading order in 1/N in this region of t. Substituting z = Swl/t — l/t c and 



using H = £wl — l/t c , we have from Eq. 



Cpf(up,u q ) — 2Nu p J dzzJo(upz)J (u q z) (62) 



= 2 N~u UpJi(u P )Jo(u q ) - UqJpjup) Ji(u q ) ,g 3 , 

where we have defined u p = Eu p and u q = Eu q and used Eq. (11.3.29) of 
Ref. []5S|. In particular, we obtain for k = 1,2 



Ps[U 



) = UJ 2 o(u) + J1(u)) (64) 
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Ps{Ul,U 2 )=U 1 U2 



VJiO^JoQa) - u 2 Jo{u l )Ji{u 2 ) 
uj - u\ 



(65) 



where we have subtracted the disconnected part of the two-point function in 
fl5"5]). Eq. (6^) is in agreement with the result of Ref. [IS]. The zero-temperature 
limit of Eq. fl6"5|) agrees with the result obtained in Ref. M. 



5.3 General case 



We now generalize the results of the previous subsection to the case where Y = 
diag(yi, . . . , 2/ at). We only require that the y n be real. The chiral condensate for 



this case has been computed in Ref. ||15|| , and the condition for the existence 
of a condensate was found to be 

1 N 1 

iv n =l tin 

Note that there is always a condensate if some of the y n are zero (or more 
precisely, if a finite fraction of the y n is zero). We shall restrict ourselves to 
choices of Y for which the above condition is satisfied. 



We are now dealing directly with the two integrals in Eqs. fl56| ) and (157]). In 
Appendix these two integrals are evaluated in saddle-point approximation 
in the large- A" limit. Again, /i (t ) /2 (*) is of order 1/N for t < t c and suppressed 
exponentially for t > t c , respectively, where now 

tc= ^hA- (67) 

Note that t c > 1 by ( |66|) and that t c is infinite if some of the y n are zero. 
Again, only the interval 1 < t < t c contributes to leading order in 1/N, and 
we obtain from Eqs. ( |B.8| ) and ( |B.15| ) using Eqs. ( |C.8| ) and 



fi(t)f 2 (t) = ^J (i:\x)J (Z\x)-^ (68) 

to leading order in 1/N in this region of t, where C 2 {x) is given by Eq. (|C.6|) 
and x is implicitly given by Eq. ( [TOD below. To proceed, we substitute z = Y?x 
and obtain 

dz = Y?^dt. (69) 
dt y ' 
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To compute dx/dt, we use the saddle-point equation ( |C2 



) in the form 



E2f ^S^ (70) 

and differentiate with respect to t to obtain 

1 N 1 dx 1 dx 

S ' = ^N^iyi + ^yTt = 'Yx^Tt ' (71) 

where we have used Eq. ( p.6|) in the last equality. Simple algebra yields 
x 2 1 

dt = -T^zdz (72) 



C'Ux) 2S 



so that Eq. flog ) in connection with Eq. flog) reduces to Eq. ( pz|) with H = 
£ 2 x(t = 1), where x(t = 1) is the solution of (70) for t = 1. We have thus 



reduced the general case to the special case considered in the previous sub- 
section. The only complication in the general case is that one has to solve 
the higher-order polynomial equation ( |T0"D for x(t = 1). This is no problem 
numerically since we have shown in Appendix |C| that there is only one real 
and positive solution which has to be determined. 



6 Discussion 



We have considered a random matrix model appropriate for the Dirac op- 
erator of QCD at finite temperature and derived the spectral density, the 
universal spectral correlations on the scale of the mean level spacing, and the 
microscopic limits of the spectral density and spectral correlation functions. 
The calculation was done for an arbitrary real diagonal matrix added to the 
random matrix of the Dirac operator. The fact that the functional form of 
the microscopic correlations remains unchanged even with an arbitrary de- 
terministic matrix added adds further evidence to the conjecture that these 
correlations (and in particular the microscopic spectral density) are universal. 

It is a well known feature of standard random matrix theory that many of 
the conceptual and formal aspects are common to all three ensembles. In 
particular, the calculation of fluctuations shows these common features since 
the scale is always set by the mean level spacing. Therefore, we conjecture 
that the universal scaling behavior which we derived for all correlations in 
the chGUE carries over to the chGOE and the chGSE. More precisely, we 
expect that, after proper rescaling with the chiral condensate, the microscopic 
spectral correlations are given by the known zero-temperature results for the 
chGOE and the chGSE, respectively 
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There are two important directions in which one would like to generalize the 
results of this paper. First, one would like to extend our results to a non- 
zero number of flavors and non-zero topological charge. This presents certain 
technical challenges to the graded eigenvalue method, but we are convinced 
that an application of the Itzykson-Zuber integral for rectangular supermatri- 
ces [32] will solve this problem. The result, however, can be conjectured now: 
After proper rescaling with the chiral condensate, we will obtain the same de- 
pendence on Nf and v as in the zero-temperature result computed in Ref. ||. 
Second, one would like to consider complex diagonal matrices instead of just 
real ones in order to model the effect of the chemical potential on the spectrum 
of the Dirac operator. For the latter case, it is important to go to Nf 7^ 0, i.e., 
to the unquenched case, as was argued by Stephanov JLTJ. At present, there 
is quite some effort in this direction []39| . 



Note: After completion of the present work, we learned that a similar calcula- 
tion was performed simultaneously and independently by Jackson, §ener, and 
Verbaarschot. 
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A Uniform convergence 



The integral over if) in Eq. fl25| ) does not converge uniformly since the matrix 
a is arbitrary complex rather than hermitian so that its eigenvalues are com- 
plex rather than real. Note that it is sufficient to consider only the commuting 
variables of a since the anticommuting variables cannot cause convergence 
problems. Hence, the a- and the ^-integration in (|25| ) cannot simply be inter- 
changed. 

The obvious way to proceed is to deform the integration contour of the a- 
variables so that uniform convergence is obtained. This was done in Ref. Jl9| 
where the one-point function was computed by an explicit re-parameterization 
of the a-variables. This entails replacing compact integration variables by 
non-compact ones. The generalization of this approach to higher correlation 
functions is not trivial. One would have to pseudo-diagonalize a using non- 
compact groups and compute the resulting integrals of the Itzykson-Zuber 
type. 
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However, there is an easier way to proceed. In the calculation of the one-point 
function, one observes that a naive interchange of the a- and -^-integrations, 
while upsetting the real part of the result, leaves the imaginary part of the 



result intact |L9j]. The imaginary part of the result is all that is needed to 
construct the one-point function. Hence, as far as the one-point function is 
concerned, we can interchange the a- and ^-integrations and still get the 
correct result. While there is no proof at present that this remarkable fea- 
ture persists for higher correlation functions, we conjecture that it does. This 
conjecture is very strongly supported by the fact that we obtain the correct 
zero-temperature result for the higher correlation functions by simply inter- 
changing the order of the integrations as we have done in the text. Since the 
convergence problem is completely unrelated to the introduction of the matrix 
Y at finite temperature, it is hard to see how the presence of Y could change 
anything in this respect. 

Of course, this problem deserves further attention. Work in this direction is 
in progress, but there is no doubt that the results presented in this paper will 
be unaffected by a clarification of this issue. 



B Saddle-point approximation for the special case 



We wish to compute the two integrals in Eqs. (^) and (|5^) in the large-iV 
limit by saddle-point approximation. An interesting phenomenon occurs since 
the integrals depend on an external parameter t. We shall see that the integrals 
can be of different order in 1/N depending on the value of t. 

Let us start with the integral in Eq. (^) since it is conceptually simpler. 
Combining terms which grow exponentially with N, we obtain 



f 2 {t) = J ds 2 s 2 J (Z 2 u q s 2 )exp(-N£ 2 (s 2 )) (B.l) 
o 

£ 2 (s 2 ) = T l 2 ts 2 2 -\og(y 2 + s 2 2 ). (B.2) 
To minimize £ 2 (s 2 ), we differentiate 

^ = 2HHs 2 - = . (B.3) 

ds 2 y 2 + si 

There are three possible solutions for s 2 , s 2 = and s 2 = ±x with 



*=WH- (a4) 
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where t c = l/(T,y) 2 . Expanding about the saddle points yields 

1 



C 2 (s 2 + 5s 2 ) = C 2 (s 2 ) + ^£' 2 \s 2 )(5s 2 ) 2 (B.5) 



with 



£ 2 (0) = - \og{y 2 ) £ 2 {±x) = 1 - t/t c + log(£ 2 t) (B.6) 

4'(0) = 2S 2 (t - t e ) C 2 (±x) = 4S 2 t(l - t/Q . (B.7) 

The sign of C 2 ' determines the location of the minimum of C 2 . For t > t c , the 
minimum is at 0. For t < t c , the minimum is at x, and x is real. The solution 
—x lies outside the interval of integration and can be discarded. We now have 
to integrate over the fluctuations. For s 2 = x, we obtain to leading order in 
1/N 



f 2 {t)=xJ Q {Y?u q x)e- NC ^ / d(5s 2 )e-^^ 




fl e - NC ^J (E 2 u q x) l X _ -. for t<t c . (B.8) 

A'O) 



For S2 = 0, we obtain to leading order in 1/N 



00 


1 -NC 2 (0) 

= ^777^- for t>t c . (B.9) 

N £ 2 '(0) v ; 

Now consider the integral in Eq. fl58|). We again combine terms which grow 
exponentially with N to obtain 



00 

hit) = lmj ds lSl I (E 2 u pSl ) exp (-JVA(si)) (B.10) 


C 1 (s l ) = T l Hsl + log(y 2 -sl), (B.ll) 

where we have omitted the imaginary increment of s± in C\. Differentiation 
yields 



^ = 2^t Sl - -p- 2 I . (B.12) 
osi y l — sf 
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This equation is identical to Eq. (p.3| ) if s\ is replaced by zs 2 . Hence, possible 
solutions are Si = and si = ±ix. Expansion about the saddle points yields 

d(0) =\og(y 2 ) C x {±ix) = t/t c - 1 - log(£ 2 t) (B.13) 

£"(0) = 2S 2 (t - t c ) £[{±ix) = 4£ 2 t(l - t/t c ) . (B.14) 

Note that £i(zs 2 ) = — £ 2 (s 2 ). Although the saddle-point equations ( |B.3| ) and 
( |B.12|) look remarkably similar, the structure of the saddle points is now quite 
different. In the following, we show that the saddle point at ix (—ix) makes 
the dominant contribution for t < t c (t > t c ). 

Consider first the case t < t c where x > 0. The possible saddle points ±ix are 
purely imaginary. Hence, we have to deform the contour of integration into 
the complex plane to reach these saddle points. Due to the positive sign of the 
imaginary increment of s\, only the point +ix can be reached. Since C"(ix) is 
positive for t < t c , we have to integrate over quadratic fluctuations parallel to 
the real axis. Since the original range of integration was to oo and since the 
saddle point lies on the imaginary axis, the integration over the fluctuations 
is only from to oo. We obtain for the contribution from the saddle point ix 



/i(i) = Im ixI (E 2 u p ix)e- NC ^ m) J d(5s 1 )e-^ c '^ Ssi)2 

o 

= iS e ' M ' (g)j °( S V) / X for t<t c . (B.15) 

V 2iV yJC'{(ix) 

The saddle point si = does not contribute at all since the imaginary part of 
the product of the value of the integrand at the saddle point and the integral 
over the fluctuations is zero in this case. 

Now consider the case t > t c where x is imaginary. The relevant saddle point 
is now Si = —ix which is real and positive. The point ix is real and negative 
and, hence, outside the interval of integration. Since C"(—ix) is negative for 
t > t c , we now have to integrate over quadratic fluctuations parallel to the 
imaginary axis from to oo. We obtain 



■JO 

/i(t) = Im {-ix)I {-J: 2 u p tx)e- NCl{ - i2) J d(i5 Sl )e-% c "(-™K i5s ^ 

o 



JL e -iV£ l( -^) Jo(S 2 -) ™ for f>tc (B16) 

2N ^J-C'K-ix) 

which is identical in structure with Eq. ( B.15|) . Again, the saddle point si = 
does not contribute for the same reason as in the case t < t c . 
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Combining Eqs. ( |B.8|) and (P.15| ) we obtain Eq. ([61] ) which is valid for t < t c . 
Hence, fi(t)f 2 (t) ~ 1/JV in this region of t. For t > t c , Eqs. flB~9D and ( |B~T6|) 
combine to yield 

mm = tm^ht/ ° SAa^a) t/te _ t (b.it) 

showing that fi(t)f2(t) is suppressed exponentially in this case. 

We would like to add that we have checked Eqs. ( |B.15|) and ( |B.16| ) as well as the 
corresponding results of Appendix |C| by an independent finite- iV calculation 
and taking the limit N — > oo. 



C Saddle-point approximation for the general case 



We now compute the more general integrals in Eqs. and ([57]) in the 

large-iV limit by saddle-point approximation. We shall encounter the same 
phenomena as in Appendix |B[ For simplicity, we restrict ourselves to the case 
where all y n are different and non-zero as in Sec J4.2| . It is straightforward to 
modify our arguments to work for special cases (some of the y n are equal or 
zero) as well. 



Again, we write Eq. fl57|) in the form (|B.2|) , where now 



1 N 

C 2 (s 2 ) = Y?ts 2 2 - _5>g(yS + si) . (C.1) 

iV n=l 

To minimize £.2(82)1 we differentiate 



dC 



2 ov2 



N 



The obvious solution S2 = persists. There are now many more solutions for 
s\ (N of them, to be precise), but we show in the following that there is only 

one real and positive solution for s\ if t < t c and no real and positive solution 

2 
2 



for s\ if t > t c , where t c is given by Eq. (B7j). Define 



N 



We seek real and positive solutions of g(x) — 0. We have 



JV 



9 ' w = 5SvSMW >0 
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lim g(x) = t > 

, < for t < t c 
g(0)=t-t c { (C.4) 

> for t > t c . 

The function g(x) increases monotonically for x > and is continuous in this 
region. It reaches a positive value as x — > oo. Thus, g(x) intersects the x-axis 
once and only once if g(0) < and not at all if g(0) > 0. This proves our 
claim. 

Let x be the only real and positive solution of g(x 2 ) = for t < t c . Thus, 
s~2 = x is the only other candidate for a minimum of Li- Expansion about the 
saddle points as in Eq. (|B.5|) yields 



1 N 1 N 

= -m E lDg ^) ^ 2 ^) = S2 ^ 2 " AT E lo S(^ + ^) ( C - 5 ) 

iV n=l iV n=l 

1 * 1 

4'(0) =2S 2 (t-t c ) ^(x)=4x 2 -g (yg + _ 2)2 . (C.6) 

The explicit form of £2(0) and £2^) will n °t be needed in the following. We 
again obtain that the minimum of £2 is at x for t < t c and at for t > t c . 
Integration over the fluctuations yields Eqs. ( p3.8|) and (|B.9| ) as before. 



Writing Eq. (^) in the form ( |B . 1 1|) yields 

1 N 

C 1 (s l ) = i:Hsl + -^og(y 2 n -sl). (C.7) 



71=1 



Note that £i(is 2 ) = —C-2(s2)- For the Si-integration, there are now many 
more candidates for saddle points. Apart from the trivial solution s\ = 
(which does not contribute for the reason discussed in Appendix 0) , there are 
2A^ other solutions for s\. Defining a function similar to g(x) and investigating 
its pole structure, we find that for t < t c there are A^ — 1 real and positive 
solutions for s\. None of these corresponds to a minimum of C\ on the real 
axis since CJ[ < at these points. The remaining two solutions are s\ = ±ix, 
and expansion about the saddle points yields 



c{(±ix) 



As in Appendix |B|, the relevant saddle point is s\ 
the fluctuations yields Eqs. ( |B.15| ) as before. Eqs. 
to yield Eq. 



ix. 



(C.8) 
(C.9) 

and integration over 



p.8| ) and ( |B.15| ) combine 
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For t > t c , there are N real solutions for (apart from si = 0). The negative 
solution for sf for t < t c has now becomes positive, and it is straightforward 
to see from the structure of the function g(x) that it is now the smallest of all 
possible solutions for s\. Note that we still have C'{ < at all possible saddle 
points. We now deform the integration contour as in Appendix [FJ: We move on 
the real axis from to the smallest positive saddle point and then up parallel 
to the imaginary axis. Only the smallest positive saddle point contributes to 
fi(t), and we obtain a result similar to Eq. ( |B.16|) , with — ix replaced by 
the smallest positive saddle point. The product /i(£)/2(i) is again suppressed 
exponentially for t > t c so that only the interval 1 < t < t c contributes to the 
integration over t in Eq. (|55|) . 
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